Fast decimeter-level GNSS positioning

ABSTRACT

Methods and apparatus are present for determining position of a rover from observations of GNSS signals. Observations of GNSS signals are obtained at a rover location, and observations of the GNSS signals are obtained at a plurality of reference stations, each reference station defining a respective baseline between the rover location and a reference station location. For each reference station, a respective differentially-corrected rover position is determined, wherein at least one of the differentially-corrected rover positions is based on one of (i) a multiple-frequency ionosphere-free observable combination, (ii) a multiple-frequency code-phase observable combination, and (iii) a single-frequency carrier-phase and code-plus-carrier-phase combination. A weighted combination of the differentially-corrected rover positions is prepared.

CROSS REFERENCE TO RELATED APPLICATIONS

This application is a continuation of U.S. patent application Ser. No.11/786,017, filed 9 Apr. 2007, now U.S. Pat. No. 7,498,979, the contentof which is incorporated herein by this reference. The content of U.S.patent application Ser. No. 10/696,528 filed Oct. 28, 2003 isincorporated herein by this reference, and the content of ProvisionalU.S. Patent Application No. 60/792,911 filed 17 Apr. 2006 isincorporated herein by this reference.

FIELD OF THE INVENTION

The present invention relates to the field of Global NavigationSatellite Systems. More particularly, the present invention relates toprocessing of signals received from multiple satellites of such systems.

BACKGROUND OF THE INVENTION

Global Navigation Satellite Systems (GNSS) include the GlobalPositioning System (GPS), the Glonass system, and the proposed Galileosystem. Each GPS satellite transmits continuously using two radiofrequencies in the L-band, referred to as L1 and L2, at respectivefrequencies of 1575.41 MHz and MHz. Two signals are transmitted on L1,one for civil users and the other for users authorized by the UnitesStates Department of Defense (DoD). One signal is transmitted on L2,intended only for DoD-authorized users. Each GPS signal has a carrier atthe L1 and L2 frequency, a pseudo-random number (PRN) code, andsatellite navigation data. Two different PRN codes are transmitted byeach satellite: a coarse acquisition (C/A) code and a precision (P/Y)code which is encrypted for DoD-authorized users. Each C/A code is aunique sequence of 1023 bits, which is repeated each millisecond.

SUMMARY OF THE INVENTION

Improved methods and apparatus for processing GNSS signals are providedby embodiments in accordance with the present invention.

PART 1: GNSS Post-processing engine. In accordance with some embodimentsof the invention, a post-processing engine is provided which enables atwo decimeter (1-sigma) horizontal positioning accuracy after anoccupation time of two minutes with dual-frequency differential carrierphase positioning for baselines up to several hundred kilometers. Thisis a great performance leap compared to the previous processing engines,where occupation time was usually several times as long for the sameaccuracy. To achieve high accuracy with low occupation time, a number ofinnovative techniques have been developed and employed, such as thesuper-smoothing of the reference station data using data from severalhours before and after the matched rover occupation time, densifying thereference observation to any interval through filtered interpolation,forward and backward filtering and final smoothing, application ofwhitening of noise technique to remove redundant multipath states toimprove observability in short time span, using the L1 carrier phase andthe L1 carrier phase plus L1 code observable combination to efficientlyfilter the single-frequency case, and detection of movement of the roverduring static periods. Accuracy can be further enhanced in accordancewith some embodiments of the invention with multi-baseline adjustment.

PART 2: Post-processing Accuracy (PPA) Predictor. In accordance withsome embodiments of the invention, post-processing accuracy is predictedduring data collection using a look-up table in which the inputvariables include one or more of: duration since carrier lock wasacquired; HDOP (horizontal dilution of precision), number of GNSSsatellites being continuously tracked; and whether single or dualfrequency carrier data is being logged. These parameters are fed into afunction that looks up a PPA value from a pre-defined table. The valuesin the table have been assigned as a result of prior processing ofsample data with said parameters. Other input variables can also beused, such as PDOP instead of HDOP, and the knowledge of whether theantenna is static or moving. Predicted accuracy in accordance with someembodiments is horizontal RMS (root mean square). The same techniquescan be used to predict vertical accuracy, for example with differentinput variables such as VDOP (vertical dilution of precision).

These and other aspects and features of the present invention will bemore readily understood from the embodiments described below withreference to the drawings, in which:

FIG. 1 is a flowchart of a post-processing engine in accordance with anembodiment of the invention;

FIG. 2 shows reference-station locations of the International GNSSService;

FIG. 3 illustrates a situation where reference data is available atepochs t₀ and t₂, and reference data at an intervening epoch t₁ isgenerated using interpolation in accordance with an embodiment of theinvention;

FIG. 4A and FIG. 4B illustrate a principle of interpolation inaccordance with some embodiments of the invention;

FIG. 5 illustrates the input and output data of an interpolationsub-module in accordance with some embodiments of the invention;

FIG. 6 is a flow chart showing a process for interpolation of referencestation data in accordance with some embodiments of the invention;

FIG. 7 illustrates the effect of a carrier-phase cycle slips causing adistinct jump in carrier-phase measurements;

FIG. 8 illustrates in more detail portions of FIG. 1 in accordance withsome embodiments of the invention;

FIG. 9 shows an exemplary embodiment of a module which forms theminimum-error-combinations for dual-frequency GNSS data in accordancewith some embodiments of the invention;

FIG. 10 is a flow chart showing operation of a Kalman filter withtoggling of position and epoch in accordance with some embodiments ofthe invention;

FIG. 11 shows an exemplary embodiment of a module which forms theminimum-error-combinations for single-frequency GNSS data in accordancewith some embodiments of the invention;

FIG. 12 shows in more detail the process of forming single differencesfor single-frequency GNSS data in accordance with some embodiments ofthe invention;

FIG. 13 illustrates the effect of forward and backward smoothing inaccordance with some embodiments of the invention;

FIG. 14 illustrates the use of smoothing for stop-and-go processing inaccordance with embodiments of the invention;

FIG. 15 illustrates the benefit of forward and backward smoothing forprocessing of kinematic data in accordance with some embodiments of theinvention;

FIG. 16 shows an overview of a multi-baseline adjustment process inaccordance with some embodiments of the invention;

FIG. 17 is an overview of single-baseline processing with automaticcode-carrier correction in accordance with some embodiments of theinvention;

FIG. 18 is an overview of multi-baseline adjustment in accordance withsome embodiments of the invention;

FIG. 19 is an overview of an exemplary process for creating a lookuptable useful for post-processing-accuracy determination in a rover inaccordance with some embodiments of the invention;

FIG. 20 shows a process for creating a PPA lookup table in accordancewith some embodiments of the invention;

FIG. 21 illustrates for each of various HDOP values (bins) aparameterized relation between RMS positioning error and duration ofcontinuously-tracked rover data (for a given number of satellitestracked) in accordance with some embodiments of the invention;

FIG. 22 illustrates for each of various numbers of satellites tracked(bins) a parameterized relation between RMS positioning error andduration continuously-tracked rover data (for a given HDOP value) inaccordance with some embodiments of the invention;

FIG. 23 shows a process using a PPA lookup table in a GNSS roverreceiver during data collection in accordance with some embodiments ofthe invention;

FIG. 24 is an operational view of GNSS data collection with a rover inaccordance with some embodiments of the invention; and

FIG. 25 is a high-level view of GNSS data-collection operations inaccordance with some embodiments of the invention

DETAILED DESCRIPTION

Part 1: GNSS Post Processing Engine

Some embodiments in accordance with the invention provide a GNSSpost-processing engine which enables a user to achieve two decimeter(1-sigma) horizontal accuracy over a short occupation time of twominutes with dual-frequency differential carrier phase positioning forbaselines up to several hundred kilometers. This is a great performanceleap compared to the previous processing engines, where occupation timewas usually several times as long for the same accuracy. To achieve thehigh accuracy, low occupation time goal, a number of innovativetechniques have been developed and employed, such as the super-smoothingof the reference station data using data from several hours before andafter the matched rover occupation time, densifying the referenceobservation to any interval through filtered interpolation, forward andbackward filtering and final smoothing, application of whitening ofnoise technique to remove redundant multipath states to improveobservability in short time span, using the L1 carrier phase and the L1carrier phase plus L1 code observable combination withsingle-differencing to efficiently filter the single-frequency case, anddetection of movement of the rover during static periods. All thesetechniques will be discussed in detail below.

Post Processing Engine Flowchart. FIG. 1 is a flowchart of apost-processing engine 100 in accordance with an embodiment of theinvention. It contains seven basic modules for post-processing a set orrover data 105 with reference data 110: (1) module 115 gets thereference data, (2) module 120 interpolates and smooths the referencedata from 30-second intervals to 1-second intervals, (3) module 125matches the reference and rover data, (4) module 130 forms thesingle-difference minimum-error-combination observations, (5) module 135prepares a backward estimate of the position and ambiguity, (6) module140 prepares a forward estimate of the position and ambiguity, and (7)module 145 smooths the positions by combining the backward and forwardestimate results. Each module is explained below.

Module 1: Get Reference Data from a Reference Station.

The field hardware is a single GPS rover, such as a Trimble GeoExplorer2005 series or Trimble GPS Pathfinder ProXH or Trimble GeoXH handheldGPS receiver, with or without optional external Zephyr antenna. The userdoes not need a reference receiver. Raw GPS measurement and navigationdata can be downloaded free of charge from many sources, such as theInternational GNSS Service (IGS) website www.igscb.jpl.nasa.gov, whichincludes more than 300 stations worldwide as shown in map 200 of FIG. 2.

Module 2: Interpolate and Smooth the Reference Data.

The raw reference-station data can be used in post-processing incombination with the rover data. However, there is a drawback of thesefree data. They are generally decimated at 30-second intervals fordata-storage purposes. Many rover receivers operate at 1-secondintervals and there is thus a data synchronization issue between thereference data and the rover data (if only matching epochs are used,much of the rover data information is discarded). To overcome thisproblem, the reference data is interpolated from 30-second intervals toone-second intervals. FIG. 3 illustrates a situation where referencedata from a reference station 300 is only available at epochs t₀ and t₂(30 seconds apart), and reference data at epoch t₁ is generated usinginterpolation in accordance with an embodiment of the invention.

A known technique for interpolating GPS data is described in L. MADER etal., “Using Interpolation and Extrapolation Techniques to Yield HighData Rates and Ionosphere Delay Estimates from Continuously OperatingGPS Networks”, ION GPS 2002, 24-27 Sep. 2002, Portland, Oreg., NationalGeodetic Survey, NOS/NOAA, Silver Spring, Md. This technique uses ahigh-order polynomial (typically 8^(th) order). This approach suffersfrom the fact that a missing reference epoch or a cycle slip wouldprevent interpolation being done for up to the number of polynomialorder epochs (if an 8^(th) order polynomial is used, then interpolationwill not be possible for up to 8 epochs of reference epochs).

In contrast, embodiments in accordance with some embodiments of thepresent invention employ an interpolation method in which the geometryis subtracted from the carrier phase and code measurement using theephemeris, the remaining geometry-free errors are interpolated (ormodeled with a Kalman filter), and the geometry is added back to theinterpolated geometry-free errors. This is possible because theremaining errors are fairly low frequency.

Equation (1) lists the observation equations for GNSS code andcarrier-phase observations for the i^(th) carrier frequency (for GPSi=1, 2 currently).

$\begin{matrix}{{\rho_{i} = {R + {\delta\; R} + {c\left( {{\delta\; t_{r}} - {\delta\; t_{s}}} \right)} + T + {\frac{\lambda_{i}^{2}}{\lambda_{1}^{2}}I} + ɛ_{\rho_{i}}}}{{\lambda_{i}\left( {\phi_{i\;} + N_{i}} \right)} = {R + {\delta\; R} + {c\left( {{\delta\; t_{r}} - {\delta\; t_{s}}} \right)} + T - {\frac{\lambda_{i}^{2}}{\lambda_{1}^{2}}I} + ɛ_{\phi_{i}}}}} & (1)\end{matrix}$Where:

ρ_(i): Code/pseudorange measurement

φ_(i): Carrier phase measurement

R: Geometric range between the receiver and the satellite

δR: Satellite orbit error

c: Vacuum speed of light

δt_(s): Satellite clock error (including the relativistic correction)

δt_(r): Receiver clock error

N_(i): Carrier phase initial ambiguity

λ_(i): Wavelength

I: Ionospheric error

T: Tropospheric error

ε_(ρ) _(i) ,ε_(φ) _(i) : Code and carrier multipath and measurementnoise

FIGS. 4A and 4B illustrate a principle of interpolation in accordancewith some embodiments of the invention. As shown in FIG. 4A, thecode/pseudorange measurements 405, 410, 415 of the reference station,e.g., at discrete 30-second intervals, are a composite of known, mainlygeometry-related elements (e.g., satellite position, reference stationposition, satellite and receiver clock error, modelled troposphericeffect, relativity) and unknown elements (e.g., satellite orbit/clockerror, residual tropospheric error, ionospheric error, multipath,noise). The known/geometry-related elements are subtracted, leaving aresidual error at each discrete time-point t_(j)

$\begin{matrix}{\mspace{20mu}{{{ɛ_{\rho_{i}}\left( t_{j} \right)} = {{\rho_{i}\left( t_{j} \right)} - {R\left( t_{j} \right)} - {c\left\lbrack {{\delta\;{t_{r}\left( t_{j} \right)}} - {\delta\;{t_{s}\left( t_{j} \right)}}} \right\rbrack} - {T\left( t_{j} \right)} - {\frac{\lambda_{i}^{2}}{\lambda_{1}^{2}}{I\left( t_{j} \right)}}}}{{ɛ_{\phi_{i}}\left( t_{j} \right)} = {{\lambda_{i}\left\lbrack {{\phi_{i\;}\left( t_{j} \right)} + N_{i}} \right\rbrack} - {R\left( t_{j} \right)} - {c\left\lbrack {{\delta\;{t_{r}\left( t_{j} \right)}} - {\delta\;{t_{s}\left( t_{j} \right)}}} \right\rbrack} - {T\left( t_{j} \right)} - {\frac{\lambda_{i}^{2}}{\lambda_{1}^{2}}{I\left( t_{j} \right)}}}}}} & (2)\end{matrix}$

A simple interpolation approach in accordance with an embodiment of theinvention is to determine a linear interpolation for example at 1-secondintervals between the discrete samples at 30-second intervals, e.g.,interpolate a residual 430 at time t between residual 435 at t₁ andresidual 440 at t₂ as illustrated at 440 in FIG. 4B.

With the phase change rate between two time points t₁ and t₂

$\begin{matrix}{{{\delta\varphi}_{i} = {{\lambda_{i}{\delta\phi}_{i}} = \frac{{ɛ_{\phi_{i}}\left( t_{2} \right)} - {ɛ_{\phi_{i}}\left( t_{1} \right)}}{t_{2} - t_{1}}}}\;} & (3)\end{matrix}$

the interpolated phase residual at an intermediate time point t isε_(φ) _(i) (t)=ε_(φ) _(i) (t ₁)+(t−t ₁)δφ_(i).  (4)

For the dual frequency case i=1, 2 the phase change rate can be splitinto two parts. A frequency independent (geometric) part

$\begin{matrix}{{\delta\varphi}_{G} = {{{- \frac{\lambda_{2}^{2}}{\lambda_{1}^{2} - \lambda_{2}^{2}}}{\delta\varphi}_{I}} + {\frac{\lambda_{1}^{2}}{\lambda_{1}^{2} - \lambda_{2}^{2}}{\delta\varphi}_{2}}}} & (5)\end{matrix}$

and a frequency dependent (ionospheric) part

$\begin{matrix}{{\delta\varphi}_{I} = {{- \frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}}{\left( {{\delta\varphi}_{1} - {\delta\varphi}_{2}} \right).}}} & (6)\end{matrix}$

With that the code change rate can be expressed as

$\begin{matrix}{{\delta\rho}_{i} = {{\delta\varphi}_{G} + {\frac{\lambda_{i}^{2}}{\lambda_{1}^{2}}{\delta\varphi}_{I}}}} & (7)\end{matrix}$

and the interpolated code residual at an intermediate time point t isε_(ρ) _(i) (t)=ε_(ρ) _(i) (t ₁)+(t−t ₁)δρ_(i).  (8)

Now, with the interpolated code- and phase residual, the code and thephase at the intermediate time point t is computed as

$\begin{matrix}{{{\rho_{i}(t)} = {{R(t)} + {c\left\lbrack {{\delta\;{t_{r}(t)}} - {\delta\;{t_{s}(t)}}} \right\rbrack} + {T(t)} + {\frac{\lambda_{i}^{2}}{\lambda_{1}^{2}}{I(t)}} + {ɛ_{\rho_{i}}(t)}}}{{\lambda_{i}\left\lbrack {{\phi_{i\;}(t)} + N_{i}} \right\rbrack} = {{R(t)} + {c\left\lbrack {{\delta\;{t_{r}(t)}} - {\delta\;{t_{s}(t)}}} \right\rbrack} + {T(t)} - {\frac{\lambda_{i}^{2}}{\lambda_{1}^{2}}{I(t)}} + {ɛ_{\phi_{i}}(t)}}}} & (9)\end{matrix}$

The advantage of this interpolation scheme is, that the code isinterpolated using the phase change rate. This is more precise than acode-only based interpolation of the code measurements.

Another approach in accordance with an embodiment of the invention is tointerpolate for example at 1-second intervals between the discretesamples at 30-second intervals using a linear Kalman filter. The linearfunction of the Kalman filter is modelled by its bias (offset and changerate). The Kalman filter filters out noise, possibly more effectivelythan a linear interpolation, as it is a weighted interpolation(adjustment). A minimum of two samples allows a linear interpolation,while more samples allows reduction of errors. The Kalman filter is runon the small residual to a time point of interest (e.g., a desired1-second interval point t) and the change rate is obtained from theKalman filter for use as a linear interpolation factor from areference-sample time to that point.

FIG. 5 shows the input and output data of the interpolation sub-module500 of interpolation and smoothing module 120. The input data 505 aredual-frequency pseudorange and carrier-phase measurements from thereference stations, typically at 30-second intervals. The output datafrom interpolation sub-module 510 are pseudorange and carrier-phasemeasurements at, for example, one-second intervals.

FIG. 6 is a flow chart showing a process 600 for interpolation ofreference station data in accordance with some embodiments of theinvention. Ephemeris data and reference-station location data 605 areknown. Reference-station data 610, e.g., pseudorange and dual-frequencycarrier measurements at 30-second intervals from IGS sources, isaccessed. The ephemeris data and reference station location are used at615 to determine a geometry component 620 of the reference-stationmeasurement at each reference-station-measurement time tRef. For eachreference-station-measurement time tRef of interest, the geometrycomponent is subtracted from the reference-station measurement at 625 toproduce a geometry-free residual 630. These residuals are supplied to aninterpolation process 635 (e.g., a linear interpolator, or a linearKalman-filter which produces a weighted interpolation) to preparegeometry-free residuals 640 at times tRover corresponding to the sampletimes of the rover measurements to be post-processed (e.g., at 1-secondintervals). The ephemeris data and reference station location are usedat 645 to determine a geometry component 650 of the reference-stationmeasurement at each rover-data-sample time tRover (e.g., at 1-secondintervals). For each rover-data-sample time tRover, the interpolatedgeometry-free residual and the determined geometry component are addedat 655 to produce an interpolated reference measurement 660.

Another problem of using public reference data is that the quality ofthe reference receiver varies significantly. Although it is assumed thatmost IGS/CORS reference stations use geodetic quality receivers, it isstill highly probable that the user will encounter suboptimal referencedata, for example, data having strong multipath due to bad antennalocation. To overcome this problem, a super-smoothing algorithm isemployed on the reference code measurements to reduce its multipathsubstantially. This smoothing can be performed using all the referencedata available, independent of the rover data. The rover may only havean occupation of two minutes for a specific point of interest; howeverthe reference-station data can be smoothed by using all thereference-station data for a period from several hours before untilseveral hours after the rover's occupation time since the smoothing hasnothing to do with the rover data itself, only the time matters.

Smoothing merges ‘absolute’ pseudorange capability and ‘relative’carrier phase capability. The GPS pseudorange measurement is noisy butnot ambiguous; the GPS carrier-phase measurement is precise butambiguous. Precise and unambiguous pseudorange measurement can beproduced by combining these properties. The method is to average thegeometry-free, ionosphere-free code-minus-carrier observable, whichessentially contains only multipath, noise, and ambiguity. If thetracking is continuous for a satellite for several hours, the noise andmultipath can be effectively removed, giving an accurate ambiguity. Thesmoothed pseudorange at any epoch is recovered by adding the averagedambiguity to the epoch value of the carrier phase combination used tocompute the ambiguity. Thus the precision of the smoothed reference codewill approach that of the carrier phase.

In the dual frequency case i=1, 2, using the observation equations (1),a divergence-free smoothing formula can be constructed as follows:

$\begin{matrix}{{{{\rho_{1} - {\lambda_{1}\phi_{1}} - {2\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}\left( {{\lambda_{1}\phi_{1}} - {\lambda_{2}\phi_{2}}} \right)}} = {{{{- \lambda_{1}}N_{1}} - {2\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}\left( {{\lambda_{1}N_{1}} - {\lambda_{2}N_{2}}} \right)}} = N_{s\; 1}}}{\rho_{2} - {\lambda_{2}\phi_{2}} - {2\frac{\lambda_{2}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}\left( {{\lambda_{1}\phi_{1}} - {\lambda_{2}\phi_{2}}} \right)}} = {{{{- \lambda_{2}}N_{2}} - {2\frac{\lambda_{2}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}\left( {{\lambda_{1}N_{1}} - {\lambda_{2}N_{2}}} \right)}} = N_{s\; 2}}}\mspace{20mu}{Where}} & (10) \\{\mspace{20mu}{\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}\left\lbrack {{\lambda_{1}\left( {\phi_{1} + N_{1}} \right)} - {\lambda_{2}\left( {\phi_{2} + N_{2}} \right)}} \right\rbrack}} & (11)\end{matrix}$is the ionospheric error on L1.

As long as there is no cycle slip, then the right side of Equation (10)is a constant. Code super-smoothing is based on this property; itaccumulates the left side of Equation (10), and calculates the meanvalues Ñ_(s1) and Ñ_(s2) over multiple epochs. The smoothed pseudo-rangeat epoch (t_(k)) can be calculated as

$\begin{matrix}{{{\rho_{s\; 1}\left( t_{k} \right)} = {{\lambda_{1}{\phi_{1}\left( t_{k} \right)}} + {2\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}\left( {{\lambda_{1}{\phi_{1}\left( t_{k} \right)}} - {\lambda_{2}{\phi_{2}\left( t_{k} \right)}}} \right)} + {\overset{\sim}{N}}_{s\; 1}}}{{\rho_{s\; 2}\left( t_{k} \right)} = {{\lambda_{2}{\phi_{2}\left( t_{k} \right)}} + {2\frac{\lambda_{2}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}\left( {{\lambda_{1}{\phi_{1}\left( t_{k} \right)}} - {\lambda_{2}{\phi_{2}\left( t_{k} \right)}}} \right)} + {\overset{\sim}{N}}_{s\; 2}}}} & (12)\end{matrix}$

If the satellite is tracked continuously for an extended period of timewithout any cycle slip, then the smoothed pseudo-range will be asprecise as the carrier phase observation. A huge reduction of noiselevel is achieved.

One advantage of performing the code smoothing using Equation (12) isthat the terms

$\begin{matrix}{{{\lambda_{1}{\phi_{1}\left( t_{k} \right)}} + {2\frac{\lambda_{1}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}\left( {{\lambda_{1}{\phi_{1}\left( t_{k} \right)}} - {\lambda_{2}{\phi_{2}\left( t_{k} \right)}}} \right)}}{and}} & (13) \\{{\lambda_{2}{\phi_{2}\left( t_{k} \right)}} + {2\frac{\lambda_{2}^{2}}{\lambda_{2}^{2} - \lambda_{1}^{2}}\left( {{\lambda_{1}{\phi_{1}\left( t_{k} \right)}} - {\lambda_{2}{\phi_{2}\left( t_{k} \right)}}} \right)}} & (14)\end{matrix}$for the current epoch in Equation (12) can be evaluated before computingthe smoothed ambiguities Ñ_(s1) and Ñ_(s2), since Ñ_(s1) and Ñ_(s2) needdata from future epochs. Only after a complete satellite pass has beenprocessed, it is necessary (and possible) to add back Ñ_(s1) and Ñ_(s2)as corrections.

Carrier-phase cycle slips can occur as a result of an interruptedcarrier-phase tracking loop, causing a distinct jump in the carrierphase measurement, as shown in FIG. 7. Carrier-phase cycle slip can bedetrimental to carrier phase positioning; for example, the abovementioned smoothing will suffer from cycle slips because every timethere is a cycle slip, the right side of Equation (10) will change andthus the averaging of the left side of Equation (10) will then have tobe reset, degrading the smoothing performance. If a cycle slip isencountered, the interpolation process also has to be reset andrestarted, causing a loss of reference data for up to three epochs (thenumber of epochs to converge the Kalman filtered interpolator) of rawreference data (e.g., reference data is lost for 90 seconds if thereference data is at 30-second intervals). It is thus important todetect and repair these cycle slips. There are many prior art techniquesfor cycle slip repair. A detailed discussion on cycle slip detection andrepair, ranging from carrier minus code, ionospheric residual and tripledifference, is found for example at B. HOFMANN-WELLENHOF, GPS Theory andPractice, 2d. Ed., 1992, at Chapter 9. TEQC is a quality control checkrun by CORS which will show cycle slips in the reference-station datafile. A rover receiver will also normally provide an indication of acycle slip.

FIG. 8 illustrates in more detail the interpolation and smoothing module120, including a smoothing sub-module 800 in accordance with someembodiments of the invention. Reference-station pseudorange measurements805 (e.g., at 30-second intervals) and reference-station carrier-phasemeasurements 810 (e.g., also at 30-second intervals) are combined at 815to produce a geometry-free and ionosphere-free code-minus-carrierobservable 820 for each carrier (L1 and L2). These observables, whichare mainly multipath, noise and carrier-ambiguity, are each averagedover many epochs (for example, over many hours) at 825 to average outmultipath and noise and produce a mean carrier-ambiguity for eachcarrier frequency at 830. The mean carrier-ambiguities 830 are combinedat 835 with the interpolated carrier-phase measurements 660(interpolated as discussed above, e.g., to 1-second intervals), toproduce a set 840 of interpolated and smoothed reference codes at timesmatching the rover measurements (e.g., 1-second intervals).

Module 3: Match the Reference and Rover Data.

FIG. 8 also shows module 125 which matches the rover data 105 and the(interpolated and smoothed) reference data 840 at each epoch, andcomputes the single differences 850. After this, differential processingcan start. The next several sections will discuss innovative techniquesdeveloped for this purpose.

Module 4: Form the Minimum-Error-Combination Observations.

The reference data usually contain dual-frequency observables, L1 andL2. However the rover can provide either dual frequency or singlefrequency observables, depending on the configurations and the hardware.(For example, the Trimble GeoExplorer 2005 series and Trimble GPSPathfinder ProXH and Trimble GeoXH handheld GPS receivers providesingle-frequency observables when used with the internal antenna andprovide dual-frequency observables when used with the optional externalZephyr antenna). Depending on the number of frequencies available (e.g.,single- or dual-frequency for L1/L2 GPS), different observationcombinations and parameterization schemes are formed. Regardless of thenumber of rover receiver frequencies, two types of observations are usedin the final filtering: one code combination and one carrier phasecombination.

Dual Frequency. If both L1 and L2 carrier-phase data are available, thentwo types of observation combinations are formed. The first is theso-called minimum-error code-minus-carrier geometry-free combination.The second is the minimum-error carrier-phase combination. These areexpressed, respectively, as:ρ_(cc) =aλ ₁φ₁ +bλ ₂φ₂−(cρ ₁+(1−c)ρ₂)φ_(cc) =aλ ₁φ₁ +bλ ₂φ₂  (15)

Where

a: combination coefficient for L1 carrier phase

b: combination coefficient for L2 carrier phase

c: combination coefficient for L1 pseudo-range

ρ_(cc): combined minimum error code minus carrier combinations

φ_(cc): combined minimum carrier combinations (with very long baseline,it tends to be ionosphere free combinations)

The coefficients are determined based on priori knowledge of the carrierphase noise, pseudo-range noise, baseline length, and ionosphericcontent. For very long baselines where the ionosphere becomes thedominant differential error source, the minimum-error combinations tendto be ionosphere-free combinations, on both carrier and code.

FIG. 9 shows an exemplary embodiment of a module 130 which forms theminimum-error-combinations for dual-frequency GNSS data and feeds theseto a Kalman filter 905 which performs the filtering. L1 code and L1phase data 910, and L2 code and L2 phase data 915 are supplied to asub-module 920 which creates a geometry-free minimum-error code-carriercombination 925. L1 phase data 930 and L2 phase data 935 are supplied toa sub-module 940 which creates a geometric minimum-error combination945. These combinations 925 and 945 are supplied to a static-mode Kalmanfilter 905 having states for

Rover Position x₁, y₁, z₁

Rover Clock 1 T₁ (rover clock during epoch 1)

Rover Clock 2 T₂ (rover clock during epoch 2)

L1 ambiguities N₁₋₁ . . . N_(1-N) (one per satellite for L1)

L2 ambiguities N₂₋₁ . . . N_(2-N) (one per satellite for L2)

The filter's rover clock state is toggled between T1 and T2 from oneepoch to the next. The ambiguities change only when reset due to cycleslip. The Kalman filter provides as outputs 950 the rover positionvalues and the ambiguities.

In an alternate embodiment, the Kalman filter of FIG. 9 is enhanced forkinematic mode by adding further states for rover position and roverclock so that the states are

Rover Position 1 x₁, y₁, z₁ (rover position during epoch 1)

Rover Clock 1 T₁ (rover clock during epoch 1)

Rover Position 2 x₂, y₂, z₂ (rover position during epoch 2)

Rover Clock 2 T₂ (rover clock during epoch 2)

L1 ambiguities N₁₋₁ . . . N_(1-N) (one per satellite for L1)

L2 ambiguities N₂₋₁ . . . N_(2-N) (one per satellite for L2)

The filter's rover clock state is toggled between T1 and T2 from oneepoch to the next. When processing a current epoch, the Kalman filterkeeps track of the rover position x, y, z for the epoch previouslyprocessed and determines whether the rover position has changed. Whenthe rover has moved from the previous epoch, the rover position statesare toggled between x₁, y₁, z₁ and x₂, y₂, z₂ to operate the filter in akinematic mode. When the rover has not moved from the previous epoch,the rover position states are not toggled and the filter is thusoperated in a static mode.

FIG. 10 is a flow chart showing operation of Kalman filter 1005 withtoggling of position and epoch in accordance with some embodiments ofthe invention. At 1005 an epoch flag (rover clock) is initialized to T1and rover position is initialized to x₁, y₁, z₁. At 1010 the epoch flagis checked and, if not set to T1 is toggled at 1015 to T1. If at 1010the epoch flag is set to T1 then it is toggled at 1020 to T2. At 1025 acheck is made of the Kalman filter's position states to determinewhether position has changed more than a predetermined amount. If no,then data for an epoch is retrieved at 1030. If yes, then at 1035 theposition is checked and, if not set to x₁, y₁, z₁ it is set at 1040 tox₁, y₁, z₁. If position is found at 1035 to be set to x₁, y₁, z₁ then at1045 it is toggled to “x₂, y₂, z₂”. After toggling the position at 1040or 1045, data for an epoch is retrieved at 1030. After retrieving datafor an epoch at 1030, the epoch of data is processed at 1050 in a Kalmanfilter. The rover position and ambiguities are stored at 1055. At 1060 acheck is made for a next epoch of data and, if available, the processreturns to 1010.

Single Frequency. Single-frequency data means that the L1 carrier phasewill be used as it is since no L2 carrier is available. However, for theL1 pseudo-range measurement, the (ρ₁+λ₁φ₁) combination is used, which isionosphere-free (but ambiguous), as contrasted with the classical(ρ₁−λ₁φ₁) combination, which is not ionosphere-free. (The L1 phase plusL1 code combination is found to be ionosphere-free since the ionosphereerror has the same magnitude but opposite sign on the L1 carrier and L1code observations.) The (ρ₁+λ₁φ₁) combination is known to be used in thecontext of single-point positioning (see Y. GAO et al., High PrecisionKinematic Positioning Using Single Dual-Frequency Receiver, THEINTERNATIONAL ARCHIVES OF THE PHOTOGRAMMETRY, REMOTE SENSING AND SPATIALINFORMATION SCIENCES, Vol. 34, Part XXX), and for triple-differencedpositioning without floating solutions (see B. Remondi et al., FinalReport: Investigation of Global Positioning System Single FrequencyHardware for the U.S. Environmental Protection Agency, EPA ReferenceDW13936132-01-0, April 1994, ten pages) but has not previously beenproposed for differential positioning as in embodiments of the presentinvention.

FIG. 11 shows an exemplary embodiment of a module 130 which forms theminimum-error-combinations for single-frequency GNSS data and feedsthese to a Kalman filter 1105 which performs the filtering.Single-differenced L1 code and L1 phase data 1110 (single differencesbetween rover data and interpolated and supersmoothed reference data)are supplied to a sub-module 1120 which prepares the geometric L1 codeand phase combination. Single-differenced L1 phase data 1115, which isgeometric, is supplied via an L1 phase sub-module 1125. The main errorin the single-differenced L1 code plus L1 phase combination is multipatherror, which is time-correlated, besides the a priori unknown(time-constant) ambiguity which is estimated by the Kalman filter. FIG.12 shows in more detail the process of forming single differences forsingle-frequency GNSS data. Interpolated and smoothed referencepseudorange data 1205 and interpolated reference L1 carrier-phase data1210 and the rover observables 105 are supplied to a process which formssingle differences 1220 for pseudorange and for L1 carrier-phase.

The Kalman filter for single-frequency rover data is similar to thegeometry filter described in U.S. patent application Ser. No.10/696,528, except that it has two observables per epoch. (The referencedata needs L2 for the supersmoothing process described below) Thesingle-differences for code and L1 carrier are formed using referencedata interpolated to match the rover data epochs. FIG. 11 illustratesformation of single differences for L1-only rover data.

The L1 code plus L1 phase differences and the L1 phase differences aresupplied to a Kalman filter having states for

Rover Position 1 x₁, y₁, z₁ (rover position during epoch 1)

Rover Clock 1 T₁ (rover clock during epoch 1)

Rover Position 2, x₂, y₂, z₂ (rover position during epoch 2)

Rover Clock 2 T₂ (rover clock during epoch 2)

L1 ambiguities N₁₋₁ . . . N_(1-N) (one per satellite for L1)

The filter's rover clock state is toggled between T₁ and T₂ from oneepoch to the next. When processing a current epoch, the Kalman filterkeeps track of the rover position x, y, z for the epoch previouslyprocessed and determines whether the rover position has changed. Whenthe rover has moved from the previous epoch, the rover position statesare toggled between “x₁, y₁, z₁” and “x₂, y₂, z₂” to operate the filterin a kinematic mode. When the rover has not moved from the previousepoch, the rover position states are not toggled and the filter is thusoperated in a static mode. The ambiguities change only when reset due tocycle slip. The Kalman filter provides as outputs the rover positionvalues and the ambiguities. The flow chart of FIG. 10 also applies fortoggling of position and epoch for the L1-only case; the Kalman filterstates for processing L1-only rover data differ from those forprocessing of L1/L2 dual-frequency rover data (for L1-only theambiguities are L1 and for dual-frequency the ambiguities areminimum-error L1/L2 combination; see U.S. patent application Ser. No.10/696,528 and L. SJÖBERG , The best linear combinations of L1 and L2frequency observables in the application of Transit/Doppler and GPS,MANUSCRIPTA GEODETICA 15, 1990, pp. 17-22).

Module 5: Prepare a Backward Estimate of the Position and Ambiguity;Module (6): Prepare a Forward Estimate of the Position and Ambiguity

GPS carrier and code observations include not only white noise, but alsocolored noise (non-Gaussian errors), for example, multipath. Multipatherror on the code observations can be up to tens of meters and on thecarrier phase can be up to several centimeters. The multipath errorstypically have a time constant of tens of seconds. To achieve thehighest processing accuracy, the multipath should be properly accountedfor. There are two ways to take care of this. The first is via stateaugmentation where an additional multipath state is estimated in theKalman filter; the second is the whitening of noise technique, where theobservation of the current epoch is differenced with the observation ofthe previous epoch via a correlation coefficient. The first method ismore expensive from a computational load perspective and is thus notpreferred. The second has proven its efficiency.

The basic principle of this whitening of noise technique is demonstratedas follows. Assuming the GPS carrier phase observation at epoch t₀ is L₀and at epoch t₁ is L₁, and the multipath on the carrier phase is foundto be a first order Gauss-Markov process having a time constant oft_(c), then the following differenced observation is found to be white:

$\begin{matrix}{\phi_{1}^{\prime} = {\phi_{1} - {{\mathbb{e}}^{\frac{t_{1} - t_{0}}{t_{c}}}\phi_{0}}}} & (16)\end{matrix}$By properly taking into account the stochastic nature of the multipath,the algorithm provides optimal position accuracies and reduces theKalman filter's convergence time by reducing the number of filter statesto 8 plus N (where N is the number of satellites) from the 4 plus 4Nstates of the conventional Kalman filter approach.

One implication of this whitening of noise technique is that since theobservation contains observation not only from this epoch, but also fromthe prior epoch, the estimated states also have to contain not only theposition states (x,y,z) from this epoch, but also position states fromprevious epoch. So the complete state vector is:(x ₁ ,y ₁ ,z ₁ ,t ₁ ,x ₀ ,y ₀ ,z ₀ ,t ₀ ,N ₁ ,N ₂ . . . N _(n))  (17)Where

x₁, y₁, z₁, t₁: rover position & rover clock error states of currentepoch

x₀, y₀, z₀, t₀: rover position & rover clock states of previous epoch

N₁, N₂ . . . N_(n): ambiguity states for satellite 1 to satellite n.

Module 7: Smoothing Positions

After the data have been processed in both forward and backward mode, anadditional estimation step can be taken on the forward/backwardprocessed results, known as smoothing. Smoothing is generally well knownas an estimation technique in GNSS. An introduction to smoothingtechniques is found in A. GELB, APPLIED OPTIMAL ESTIMATION, MIT Press,1974, at Chapter 5, Optimal Linear Smoothing.

Although a goal is to achieve two-decimeter accuracy with two minutes ofdata collected at a rover position, it is rare that a user will log onlytwo minutes of data during a survey. Often what happens in the field isthat the user will occupy a feature for some time, such as two minutes,and then keep the satellite tracking (which is important, as it keepsthe ambiguity information) while moving to the next feature point to belogged. As long as satellite tracking is not interrupted, there is apossibility to use the future ambiguity information (ambiguities of asucceeding epoch) for current position estimate (processing of a currentepoch) during post-processing. The smoothing idea is based on thiscontinuous satellite and ambiguity tracking. It uses all the informationavailable from the entire survey, not just the data up to the currentepoch.

FIG. 13 illustrates the effect of smoothing over a period from to tot_(e). Curve 1305 shows the convergence of the position RMS with respectto the satellite tracking time in forward filtering; the RMS decreaseswith more data (longer satellite tracking time). The result of theforward processing isC _(pf) ,C _(af) ,C _(paf) ,p _(f) ,a _(f)  (18)

where

a_(f): forward ambiguity estimate

p_(f): forward position estimate

C_(af): forward ambiguity variance-covariance (VC) matrix

C_(pf): forward position VC matrix

C_(paf) forward position ambiguity covariance matrix

Curve 1310 shows the RMS convergence when the data is processed in thereverse direction (backward filtering). The RMS also decreases with moredata. The result of the back processing isC _(pb) ,C _(ab) ,C _(pab) ,p _(b) ,a _(b)  (19)

where

a_(b): backward ambiguity estimate

p_(b): backward position estimate

C_(ab): backward ambiguity variance-covariance (VC) matrix

C_(pb): backward position VC matrix

C_(pab): backward position ambiguity covariance matrix

The almost level line 1315 shows the smoothed RMS. The smoothing is doneas follows:a _(s)=(C _(ab) ⁻¹ +C _(af) ⁻¹)(C _(af) ⁻¹ a _(b) +C _(ab) ⁻¹ a _(f))p _(s) =p _(b) −C _(pab) ′C _(ab) ⁻¹(a _(b) −a _(s))  (20)

Where

a_(s): as smoothed ambiguity estimate

p_(s): smoothed position estimate

The advantage of smoothing is demonstrated by the fact that namely thesmoothed position precisionC _(ps) <C _(pf) ,C _(pb)  (21)

Where

C_(ps): smoothed position VC matrix

is better than both forward and back processing alone. Plot 1315 showsthat the smoothed position accuracy is substantially evenly distributedfor the entire survey, because it uses all the ambiguity information(before and after the epoch currently being processed), not just pastambiguity information as in real-time processing. This offers a bigadvantage to the user if satellite tracking is not interrupted whenmoving between feature points. The user may only occupy a feature pointfor two minutes near the beginning of the survey, e.g., at a time t₁,but if the survey is continued for a substantial time (such as one hour)without interrupting satellite tracking, the computed position for thatfeature point after the smoothing will be more accurate than with justthe forward-filtering results.

The smoothing also supports stop-and-go processing. FIG. 14 illustratesthe scenario. In the time interval from point A to point B the rover isin motion, the curve in this region representing a kinematic session.From point B to point C the rover is static. From point C to point D,the rover is in motion again, representing a kinematic session. Firstthe backward processing is done. During the static occupation periodfrom C to B, only ambiguity, position information (state estimate,variance-covariance matrix) at epoch B is stored, all previous results(from C to B) are discarded. Next, the forward processing is performed,starting from epoch A to B. Before processing data at epoch B, thesmoothing can be done with the stored backward processing results(position, ambiguity, states and Variance-Covariance matrix); theforward processing results from applying Equation (20).

Detection of Rover Movement During Static Periods:

Detection of rover movement is possible during periods reported by theuser as static. In many cases the user provides wrong stop-and-goinformation. This means that the rover is moving during a period inwhich the user has reported the rover to be static. Naturally this willcause a large error in the resulting static position computed with thepost-processing engine. To solve this problem, rover movement can bedetected, by checking the given stop-and-go intervals for staticoutliers before the engine processes these intervals in static mode.

A straightforward approach would be to compute all positions inkinematic mode. Then the differences of consecutive epochs would beinspected to detect movements. As the positions very much depend on theambiguity estimates at the specific epoch (influenced by filterconvergence and rising/setting/slipping satellites), the accuracy ofthose differences is generally not very good (decimeter-range as shownin by the upper curve 1505 in FIG. 15). To become independent ofambiguity estimates, cycle slips and constellation changes, anotherprior art method is the so called “delta-phase method”. Here, theposition difference between two epochs is estimated directly from thephase measurement differences of the same two epochs. The differencecontains no ambiguity if there has been no cycle slip.

In accordance with embodiments of the invention, the Kalman filterstates already include the rover positions of the current epoch and theprevious epoch (see Equation (17) above). These are provided for the“whitening of noise” technique to model multipath and residualionosphere as discussed above. These positions can be used to directlycompute the estimated motion of the rover between the two latest epochs,both for the forward- and the backward-processing (see the two lowerplots of FIG. 10), using the relation

$\begin{matrix}{\overset{.}{r} = \frac{{{\overset{\rightarrow}{r}\left( t_{i} \right)} - {\overset{\rightarrow}{r}\left( t_{i - 1} \right)}}}{t_{i} - t_{i - 1}}} & (22)\end{matrix}$where {right arrow over (r)}(t_(i))=(x_(i),y_(i),z_(i)) is the positionat the current time t_(i) and {right arrow over(r)}=(t_(i-1),x_(i-1),y_(i-1),z_(i-1)) is the position at the previoustime t_(i-1). The current position is computed from the previousposition using the decorrelated delta-phase observableφ(t_(i))−aφ(t_(i-1)), with the decorrelation-coefficient α.Comparison with the upper plot of FIG. 15 shows the great improvement insensitivity of this motion detection method. Velocities of about 5 cm/sor greater can be detected and reliably separated from the noisybackground. If in a given stop-and-go interval {dot over (r)}>0.3 m/s(for example as in FIG. 15) the interval is marked as invalid ornon-static. The entire interval is then processed in kinematic mode, orthe longest static sub-interval is determined and processed in staticmode.

Multi Baseline Adjustment (MBA)

In accordance with some embodiments of the invention, MBA improvesdifferential-correction accuracy by averaging the position-determinationresults from several different baselines. This technique can providehigher accuracy than any single baseline. It also has several advantagesover typical network modeling methods: it works well over long baselinedistances, is more tolerant of imprecise base coordinates, and is suitedto ad-hoc combinations of base stations.

Disadvantages of VRS for long baselines (>300 km). Differentialcorrection is performed by applying corrections to the satellitemeasurements collected at the roving GNSS (e.g. GPS) receiver. Thesecorrections are determined in one of two ways. A single, stationary GPSbase station at a known coordinate can be used to calculate thecorrections. Or a network modeling method may be used, which uses anumber of base stations around the area, and uses their knowncoordinates and their satellite measurements to create a model of thevarious error components over the relevant period of time. This model isthen used to interpolate the corrections at the rover positions, whichare then applied. This network modeling technique is also known as a VRS(Virtual Reference Station).

A VRS provides more accurate corrections than a single base station.However there are some disadvantages. Usually the base stations cannotbe too far apart: 70 km is the limit for RTK-VRS (real-time kinematicvirtual reference station), and 300 km for DGPS-VRS (differential GPSvirtual reference station). The reference positions of the base stationsmust be very accurately known; otherwise the model can fail toinitialize. The rover should be well inside the coverage area of thebase station model, but this is not always possible, for example on acoast. Also, the network of base stations takes time to configure andinitialize, and so is not suited to ad-hoc use. Lastly, the modelingrequires preferably 5 base stations, but in many locations there arefewer available. Multi-baseline adjustment has none of thosedisadvantages.

MBA method. Given m reference stations, for each reference station i={1,. . . , m} a corrected position with the coordinates (φ_(i), λ_(i),h_(i)) can be computed using an embodiment of the post-processing enginedescribed herein. To reduce the deviation from the true position (φ₀,λ₀, h₀), these corrected positions can be averaged, using a weightingfunction proportional to the quality of the positions. Usually thequality depends on the number of satellites tracked (n_(sat)), PDOP(positional dilution of precision), distance (baseline-length), numberof epochs of valid data, RMS of the residuals, etc. Some or all of thesequantities are used in the post-processing engine to estimate thecovariance for each coordinate (see Equation (21)). A good choice forthe weighting function is thus the inverse covariance 1/σ². The averagedlatitude φ(t_(j)), longitude λ(t_(j)) and height h(t_(j)) at the epocht_(j) are then

$\begin{matrix}{{{\overset{\_}{\varphi}\left( t_{j} \right)} = {\frac{1}{\sum\limits_{i = 1}^{m}w_{i}}{\sum\limits_{i = 1}^{m}{w_{i}{\varphi_{i}\left( t_{j} \right)}}}}},{w_{i} = \frac{1}{{\sigma_{i}\left( \varphi_{i} \right)}^{2}}}} & (23) \\{{{\overset{\_}{\lambda}\left( t_{j} \right)} = {\frac{1}{\sum\limits_{i = 1}^{m}w_{i}}{\sum\limits_{i = 1}^{m}{w_{i}{\lambda_{i}\left( t_{j} \right)}}}}},{w_{i} = \frac{1}{{\sigma_{i}\left( \lambda_{i} \right)}^{2}}}} & (24) \\{{{\overset{\_}{h}\left( t_{j} \right)} = {\frac{1}{\sum\limits_{i = 1}^{m}w_{i}}{\sum\limits_{i = 1}^{m}{w_{i}{h_{i}\left( t_{j} \right)}}}}},{w_{i} = \frac{1}{{\sigma_{i}\left( h_{i} \right)}^{2}}}} & (25)\end{matrix}$

FIG. 16 shows an overview of an MBA process 1600 in accordance with someembodiments of the invention. At 1605 reference data 1610 from multiplereference stations is retrieved for post-processing with rover data1615. At 1620 a single-baseline process is repeated for each of thereference stations to prepare a set of smoothed positions 1625 for eachreference station. Smoothed positions 1625 are then averaged in an MBAprocess 1630 to produce adjusted positions 1645.

FIG. 17 is an overview of single-baseline processing with automaticcode-carrier correction 1700 in accordance with some embodiments of theinvention. Reference-station data 1705 and rover data 1710 (L1 oneand/or L1/L2) are supplied to a code processor which produces anintermediate code solution 1720. Reference-station data 1705 and roverdata 1710 (L1 one and/or L1/L2) are supplied to a float-solutionprocessor 1720 which produces an intermediate L1 float solution 1725and, if dual-frequency rover data is available, an intermediate L1/L2float solution 1730. A module 1735 combines the results of 1725 and 1730by selecting the best solution at each epoch to produce an intermediatefloat solution 1740.

FIG. 18 is an overview of multi-baseline adjustment 1800 in accordancewith some embodiments of the invention. Data 1805 from a first referencestation and rover data 1810 (L1 one and/or L1/L2) are supplied to aprocessor 1815 which produces an intermediate code solution 1820 forreference station 1 and an intermediate float solution 1825 forreference station 1. Data 1830 from a second reference station and roverdata 1810 (L1 one and/or L1/L2) are supplied to a processor 1835 whichproduces an intermediate code solution 1840 for reference station 2 andan intermediate float solution 1845 for reference station 2. Data 1850from an X-th reference station and rover data 1810 (L1 one and/or L1/L2)are supplied to a processor 1855 which produces an intermediate codesolution 1860 for reference station X and an intermediate float solution1865 for reference station X. Intermediate code solutions 1820, 1840 and1860 are supplied to an MBA process 1870 to produce an adjusted codesolution 1875. Intermediate float solutions 1825, 1845 and 1865 aresupplied to an MBA process 1880 to produce an adjusted float solution1885. A module 1890 combines the best solutions from adjusted codesolution 1875 and adjusted float solution 1885 by selecting the bestsolution at each epoch to produce a final solution 1895.

Part 2: Post-Processing Accuracy Predictor (PPA)

When collecting GNSS data, it is often useful to know the accuracy ofthe determined position. Existing GPS receivers output an estimate ofthe accuracy of the position calculated at each epoch (the ‘real-timeaccuracy estimate’). But many users perform post-processed differentialcorrection to improve the accuracy of the GPS positions. Therefore, thereal-time accuracy estimate is not a good indication of the accuracythat will be achieved later. For example, the receiver may becalculating autonomous positions, or may be applying differentialcorrections to code solutions, whereas the post-processor may becalculating a more accurate carrier-phase solution.

Some embodiments of the invention solve this problem by supplying, attime of measurement, a prediction of the accuracy that will be achievedlater after post-processing. This is very useful to the user to optimizeproductivity when a certain threshold of accuracy is required. It isvery costly to re-visit a site if required accuracy is not achieved;likewise it is unproductive to collect GPS measurements at a point forlonger than necessary.

Post-processing accuracy (PPA) prediction in accordance with embodimentsof the invention is particularly useful when carrier float solutions areto be calculated during post-processing. Their accuracy increasesaccording to the duration for which measurements are continuouslycollected (in contrast to code solutions which have roughly similaraccuracy regardless of tracking duration). For example, the PredictedPost-processed Accuracy indicator answers the question “How long do Ihave to stand here to get the 25 cm post-processed accuracy that my bossrequires?”.

In accordance with some embodiments of the invention, post-processingaccuracy is predicted during data collection using a look-up table inwhich the input variables include one or more of: duration since carrierlock was acquired; HDOP (horizontal dilution of precision), number ofGNSS satellites being continuously tracked; and whether single or dualfrequency carrier data is being logged. These parameters are fed into afunction that looks up a PPA value from a pre-defined table. The valuesin the table have been assigned as a result of prior test processingwith said parameters. Other input variables can also be used, such asPDOP instead of HDOP, and the knowledge of whether the antenna is staticor moving. Predicted accuracy in accordance with some embodiments ishorizontal RMS (root mean square). The same techniques can be used topredict vertical accuracy, for example with different input variablessuch as VDOP (vertical dilution of precision).

The above methods are primarily aimed at predicting post-processed floatsolutions, which increase in accuracy according to collection duration.But the concept can be generalized to work with other types ofdifferential correction, such as code solutions. Code solutions tend tobe independent of the duration of measurement. So the lookup tablevalues will vary less over time, or may be simplified to omit theduration variable altogether.

The above methods predict the accuracy of a single baseline, i.e., usinga single GNSS base station for differential correction. But accuracy canbe improved by using data from multiple base stations in an averagingtechnique or an area network model. Thus the accuracy prediction can befurther improved by adding additional inputs about the post-processingthat will be performed later, such as whether single or multiplebaselines will be used, their number, their distance from the rover GNSSreceiver, and their geometry around the rover GNSS receiver.

The receiver calculates during data collection an estimate of theaccuracy likely to be achieved after post-processing. This helps tooptimize productivity when collecting GNSS data for which post-processedaccuracy is important. For example, the predictor answers the question“How long do I have to stand here to get the 20 cm post-processedaccuracy I need?” The predictor examines the quality of carriermeasurements and estimates how well the post-processed float solutionwill converge in the time since carrier lock was obtained.

In one embodiment, a PPA predictor program uses a lookup table approach.GNSS data is processed using the post-processing engine, and thepositioning errors are binned according to feature occupation time, HDOPand # of satellites, separately for the single- and the dual-frequencycase. For each bin, a least squares fit is then performed to generate ananalytical function that best fits the predicted position error, thefeature occupation time, HDOP and the number of satellites tracked,separately for the single- and the dual-frequency case. A PPA predictorprogram having a lookup table based on this analytical function runs,for example, in the receiver data controller of a rover receiver such asthe Trimble GeoExplorer 2005 series handheld GPS receiver or otherreceivers mentioned above.

FIG. 19 is an overview of an exemplary process 1900 for creating alookup table useful for PPA determination in a rover. GNSS measurementdata 1905 is collected in advance with the rover at a known position foran extended period (e.g., 24 hours) and this data is stored for use inpreparing the lookup table. The data preferably includes L1 and L2observations using, e.g., a Trimble handheld receiver with externalantenna. The lookup table is prepared by selecting at 1910 a subset ofthe data, using this subset to compute at 1915 differentially-correctedposition and position error and HRMS (e.g., using reference-station dataand a PP engine as described above), and binning the results at 1920.Steps 1910, 1915 and 1920 are repeated for multiple subsets of the datato produced a set of binned results 1925. The binned results areparameterized at 1930 by computing an RMS (or HRMS) error for each bin.The parameterized results are stored in a PPA lookup table 1935.

FIG. 20 shows a more detailed view of a process 2000 for creating a PPAlookup table. At 2005, the number of frequencies is selected, e.g., L1single frequency or L1/L2 dual frequency. At 2010, the time-step size isselected, e.g., the number of epochs of data to be used in computingposition. At 2015 a start time is selected for the data to be used incomputing position, e.g., the epoch of data from which to start theposition computation. At 2020, position is then computed from a subsetof rover data 2025 (e.g. 24 hours of previously-collected rover data),the subset being defined by the number of frequencies, the number ofepochs, and the starting epoch. At 2030, error between computed positionand known position is then computed, e.g., HRMS (horizontal root meansquare) error is computed from the latitude and longitude errors.(Alternately, height may also be included in computing athree-dimensional RMS position error if desired.) At 2035 the resultingrecord is binned according to number of frequencies, number of epochs,number of satellites observed during the processed interval, and HDOP(or PDOP if three-dimensional). Bin sizes are, for example, 0.5 for HDOPand 25 seconds for time (e.g., 25 epochs of 1 second per epoch).

The process is repeated at 2040 for multiple starting times (startingepoch) so that records are obtained for multiple numbers of satellitestracked and various HDOP values. At 2045 the step size (number ofepochs) is incremented, and the processing is repeated for multiplestarting times. Step sizes can vary, for example, from 10 seconds (10epochs) to 600 seconds (600 epochs). When all records have been acquiredfor one frequency setting (e.g. L1 single frequency), then the processrepeats for another frequency setting (e.g., L1/L2 dual frequency). Itwill be seen that this process generates a large set of binned records2055. The following table illustrates an exemplary structure of thebinned records:

HDOP i-th Lat Lon Height Lat Error Lon Error Height Error Freq EpochsSats (binned) Result φ₁ λ₁ h₁ φ₁-φ₀ λ₁-λ₀ h₁-h₀ HRMS 1 120 7 1.2 1 48.112.3 600.5 . . . . . . . . . → 20 cm . . . 120 7 1.4 2 48.2 12.1 600.4 .. . . . . . . . . . . 120 7 1.4 3 48.1 12.2 600.3 . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1207 1.3 n 48.3 12.1 600.4 . . . . . . . . . 1 120 6 1.8 1 48.1 12.3 600.5. . . . . . . . . → 25 cm . . . 120 6 1.7 2 48.2 12.1 600.4 . . . . . .. . . . . . 120 6 1.9 3 48.1 12.2 600.3 . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6 1.8 n48.3 12.1 600.4 . . . . . . . . . 1 120 5 2.1 1 48.1 12.3 600.5 . . . .. . . . . → 30 cm . . . 120 5 2.4 2 48.2 12.1 600.4 . . . . . . . . . .. . 120 5 2.3 3 48.1 12.2 600.3 . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . 120 5 2.2 n 48.3 12.1600.4 . . . . . . . . .

Each bin contains a lot of single position fixes (records) of thepreviously-collected rover data. For these binned positions the HRMS(horizontal error) is determined as

$\begin{matrix}{{H\; R\; M\; S} = \sqrt{\frac{1}{n - 1}{\sum\limits_{i = 1}^{n}\left\lbrack {\left( {\varphi_{i} - \varphi_{0}} \right)^{2} + \left( {\lambda_{i} - \lambda_{0}} \right)^{2}} \right\rbrack}}} & (26)\end{matrix}$where n is the number of records in a bin, φ_(i) and λ_(i) are thelatitude and longitude of the position (computed from the data subset)of the i^(th)-record of a bin, φ₀ and λ₀ are latitude and longitude ofthe known rover position.

The HRMS results of each bin are then parameterized, e.g. by fitting acurve to describe the relationship between HRMS and time (number ofepochs). For example, for a given number of satellites tracked and HDOPvalue, parametric values a₀, a₁, a₂ are determined to relate HRMS totime t (number of epochs) according to the relation:

$\begin{matrix}{{H\; R\; M\; S} = {a_{0} + \frac{a_{1}}{1 + {a_{s}t}}}} & (27)\end{matrix}$The parameters a₀, a₁, a₂ are thus related to the number of frequencies,number of satellites tracked and HDOP. The fitting function is repeatedfor each combination (of number of frequencies, number of satellitestracked, and HDOP) to obtain a set of parameters a₀, a₁, a₂ describingthe time dependence of HDOP. Parameters a₀, a₁, a₂ are place in a lookuptable which is stored in the rover for PPA determination during datacollection. FIG. 21 illustrates for each of various HDOP values (bins) aparameterized relation between RMS positioning error and durationcontinuously-tracked rover data (for a given number of satellitestracked). FIG. 22 illustrates for each of various numbers of satellitestracked (bins) a parameterized relation between RMS positioning errorand duration continuously-tracked rover data (for a given HDOP value).The relationships for HDOP values and number of satellites tracked areshown separately in FIGS. 21 and 22 for convenience of illustration. Itwill be understood that for purposes of the PPA lookup table aparameterized relation can be established for each combination of HDOPvalue and number of satellites tracked and number of frequencies.

FIG. 23 shows a process 2300 using a PPA lookup table in a GNSS roverreceiver during data collection in the field in accordance with someembodiments of the invention. When the receiver starts operation, adetermination is made at 2305 as to the number of frequencies beingtracked, e.g., whether L1 single-frequency or L1/L2 dual frequency isdetermined in Trimble handheld GNSS receivers according to the type ofantenna in use. GNSS data is acquired at 2310 and, typically, wheneveran epoch of data is available a position is computed at 2320 anddisplayed at 2325 on a display screen 2330 or the like. Processing inthe rover can include a coarse (not differentially-corrected) positiondetermination and display of position. In accordance with someembodiments of the invention, the collected data is examined at 2335(e.g., at each epoch) for the number of satellites tracked, HDOP andnumber of epochs of data continously collected. For these values, a setof parameters a₀, a₁, a₂ is looked up in the lookup table, and thecorresponding HRMS value is computed using Equation (27). This value isdisplayed as the expected PPA on a display screen 2360 or the like.Alternatively, or in addition to display of the HRMS value, a time(number of epochs) is calculated from the inverted Equation (27) for adesired HRMS value. The remaining time to achieve the desired HRMS valueis displayed so that the user knows when enough data has been collectedand/or another form of indication (e.g., an audible or visible signal)is given to indicate when enough data has been collected to achieve thedesired HRMS value.

FIG. 24 an operational view of GNSS data collection with a rover inaccordance with some embodiments of the invention. A user positions therover antenna at a feature and initiates data collection at the feature.The rover collects GNSS data and stores it to memory. The collected GNSSdata is used to look up a PPA value in a PPA lookup table, and the PPAvalue (and/or time remaining to a desired PPA value) is displayed. GNSSdata collection continues until the PPA value meets requirements. Whenthe PPA value meets requirements, the user is so informed by anend-feature display or other indication. The user terminates operationat the feature and repositions the antenna at another feature.

FIG. 25 is a further high-level view of operations in accordance withsome embodiments of the invention. A rover 2510 is used to collect datafrom satellites SV1, SV2, . . . SVN. The data from rover 2510 isoptionally transferred to a field computer 2515 such as a Pocket PC orlaptop personal computer, which may include a field data storage device2520 such as a magnetic or optical disk drive or non-volatile memory orother type of storage. The collected data is transferred by any suitablemeans to a post-processing computer 2525. Post-processing computer 2525accesses reference data from reference-data storage 2530 (e.g., local ornetworked storage). The reference data is provided by one or morereference receivers 2535. Post-processing computer 2525 preferablyemploys an embodiment of the post-processing engine described above,with or without multi-baseline adjustment as described above. Rover 2510can be used to collect a set of data from which post-processing computer2525 prepares a PPA table as described above. Such a PPA table can beused in rover 2510 or field computer 2515 during data collection todetermine a PPA value or a remaining time to obtain a desired PPaccuracy as described above.

Those of ordinary skill in the art will realize that the detaileddescription of embodiments of the present invention is illustrative onlyand is not intended to be in any way limiting. Other embodiments of thepresent invention will readily suggest themselves to such skilledpersons having the benefit of this disclosure. For example, while aminimum-error combination is employed in the examples, those of skill inthe art will recognized that many combinations are possible and that acombination other than a minimum-error combination can produceacceptable if less than optimum results; thus the claims are notintended to be limited to minimum-error combinations other than whereexpressly called for. Reference is made in detail to implementations ofthe present invention as illustrated in the accompanying drawings. Thesame reference indicators are used throughout the drawings and thefollowing detailed description to refer to the same or like parts.

In the interest of clarity, not all of the routine features of theimplementations described herein are shown and described. It will beappreciated that in the development of any such actual implementation,numerous implementation-specific decisions must be made to achieve thedeveloper's specific goals, such as compliance with application- andbusiness-related constraints, and that these specific goals will varyfrom one implementation to another and from one developer to another.Moreover, it will be appreciated that such a development effort might becomplex and time-consuming, but would nevertheless be a routineundertaking of engineering for those of ordinary skill in the art havingthe benefit of this disclosure.

In accordance with embodiments of the present invention, the components,process steps and/or data structures may be implemented using varioustypes of operating systems (OS), computer platforms, firmware, computerprograms, computer languages and/or general-purpose machines. Themethods can be run as a programmed process running on processingcircuitry. The processing circuitry can take the form of numerouscombinations of processors and operating systems, or a stand-alonedevice. The processes can be implemented as instructions executed bysuch hardware, by hardware alone, or by any combination thereof. Thesoftware may be stored on a program storage device readable by amachine. Computational elements, such as filters and banks of filters,can be readily implemented using an object-oriented programming languagesuch that each required filter is instantiated as needed.

Those of skill in the art will recognize that devices of a lessgeneral-purpose nature, such as hardwired devices, field programmablelogic devices (FPLDs), including field programmable gate arrays (FPGAs)and complex programmable logic devices (CPLDs), application specificintegrated circuits (ASICs), or the like, may also be used withoutdeparting from the scope and spirit of the inventive concepts disclosedherein.

In accordance with an embodiment of the present invention, the methodsmay be implemented on a data processing computer such as a personalcomputer, workstation computer, mainframe computer, or high-performanceserver running an OS such as Microsoft® Windows® XP and Windows® 2000,available from Microsoft Corporation of Redmond, Wash., or Solaris®available from Sun Microsystems, Inc. of Santa Clara, Calif., or variousversions of the Unix operating system such as Linux available from anumber of vendors. The methods may also be implemented on amultiple-processor system, or in a computing environment includingvarious peripherals such as input devices, output devices, displays,pointing devices, memories, storage devices, media interfaces fortransferring data to and from the processor(s), and the like. Such acomputer system or computing environment may be networked locally, orover the Internet.

Additional References incorporated herein by this reference:

-   1. K. CHEN et al., Real-Time Precise Point Positioning Using Single    Frequency Data, ION GNSS 18^(TH) INTERNATIONAL TECHNICAL MEETING OF    THE SATELLITE DIVISION, 13-16 Sep. 2005, Long Beach, Calif., pp.    1514-1523.-   2. E. GILL et al., High-Precision Onboard Orbit Determination for    Small Satellites—the GPS-Based XNS on X-SAT, 6^(TH) SYMPOSIUM ON    SMALL SATELLITES SYSTEMS AND SERVICES, September 20-24, La Rochelle,    France, 2004, pp. 1-6.-   3. S. LEUNG et al., High Precision Real-Time Navigation for    Spacecraft Formation Flying, ION GPS/GNSS 2003, 9-12 Sep. 2003,    Portland, Oreg., pp. 2182-2193.-   4. B. REMONDI et al., F INAL REPORT: INVESTIGATION OF GLOBAL    POSITIONING SYSTEM SINGLE FREQUENCY HARDWARE FOR THE U.S.    ENVIRONMENTAL PROTECTION AGENCY, EPA Reference DW13936132-01-0,    April 1994, ten pages.-   5. T. YUNCK, Single-Frequency Precise Orbit Determination, GLOBAL    POSITIONING SYSTEM: THEORY AND APPLICATIONS VOLUME II, Volume 164    Progress in Astronautics and Aeronautics, pp. 581-584.-   6. GETTING STARTED GUIDE: GEOEXPLORER 2005 SERIES, Version 1.00,    Revision A, Part Number 46506-40-ENG, October 2005, Trimble    Navigation Limited, 114 pages.-   7. WHITE PAPER: H-STAR TECHNOLOGY EXPLAINED, Trimble Navigation    Limited, 2005, 9 pages.-   8. QUICK REFERENCE GUIDE: TRIMBLE GPSCORRECT EXTENSION, Version    1.10, Revision A, April 2005, Trimble Navigation Limited, 114 pages.-   9. TRIMBLE GPSCORRECT EXTENSION: RELEASE NOTES, Version 1.11,    Revision A, October 2005, Trimble Navigation Limited, 38 pages.-   10. GETTING STARTED GUIDE: GPS PATHFINDER OFFICE SOFTWARE, Version    3.10, Revision A, Part Number 34231-31-ENG, April 2005, Trimble    Navigation Limited, 158 pages.-   11. GPS PATHFINDER OFFICE RELEASE NOTES, Version 3.10, Part Number    34195-31-ENG, Revision A, April 2005, Trimble Navigation Limited, 24    pages.-   12. USER GUIDE: GPS PATHFINDER PRO SERIES, Version 1.00, Revision A,    May 2005, Trimble Navigation Limited, 46 pages.-   13. GPS PATHFINDER TOOLS SDK: RELEASE NOTES, Version 2.00, Part    Number 37237-61-ENG, Revision A, June 2005, 44 pages.-   14. REFERENCE MANUAL: TERRASYNC SOFTWARE, Version 2.50, Revision A,    April 2005, Trimble Navigation Limited, 352 pages.-   15. GETTING STARTED GUIDE: TRIMBLE GPS ANALYST EXTENSION, Version    1.10, Revision A, Part Number 52652-01, April 2005, 240 pages.-   TRIMBLE GPS ANALYST EXTENSION: RELEASE NOTES, Version 1.10, Part    Number 52652-01, Revision A, April 2005, 28 pages.    Innovative concepts described and illustrated herein include the    following:-   1. A method of preparing GNSS reference data for processing data    acquired by a GNSS receiver at times corresponding to first    intervals, comprising:    -   a. Accessing reference data having measurements taken at times        corresponding to second intervals longer than the first        intervals,    -   b. Determining a geometric component for each measurement,    -   c. Producing for each measurement a geometry-free residual by        subtracting the geometric component from the measurement,    -   d. Interpolating the geometry-free residuals to produce        interpolated residuals for times corresponding to the first        intervals,    -   e. Determining a geometric component for each time corresponding        to an interpolated residual,    -   f. For each time corresponding to the first intervals, combining        the interpolated residual for that time with the geometric        component for that time to produce an interpolated reference        measurement for that time.-   2. The method of 1, wherein interpolating the geometry-free residual    for given time comprises preparing a linear interpolation over a    second interval spanning the given time.-   3. The method of 2, wherein preparing a linear interpolation    comprises determining phase change rate of the residuals.-   4. The method of 2, wherein the measurements are taken from GNSS    signals of at least two carrier frequencies, and wherein preparing a    linear interpolation comprises determining a phase-change rate and    determining an interpolated code residual using the phase-change    rate.-   5. The method of 1, wherein interpolating the geometry-free residual    for a given time comprises running a linear Kalman filter over    period including the given time.-   6. The method of one of 1-5, wherein the reference data comprises    carrier-phase measurements, and wherein the method comprises    performing steps 1.a.-1.e. for the carrier-phase measurements.-   7. The method of one of 1-6, wherein the reference data comprises    code measurements, and wherein the method comprises performing steps    1.a.-1.e. for the code measurements.-   8. A method of preparing GNSS reference data for processing data    acquired by a GNSS receiver during an acquisition period,    comprising,    -   a. Accessing GNSS reference data acquired during a reference        period substantially longer than the acquisition period, the        GNSS reference data comprising code measurements and        carrier-phase measurements,    -   b. Preparing a geometry-free and ionosphere-free        code-minus-carrier observable from GNSS reference data,    -   c. Averaging the code-minus-carrier observable over the        reference period to obtain respective mean carrier-ambiguity        values for each of a plurality of carrier frequencies,    -   d. Producing a set of pre-processed reference data by combining        the mean carrier-ambiguity values with the carrier-phase        measurements.-   9. The method of 8, further comprising interpolating carrier-phase    measurements as in one of 1-7, wherein producing a set of    pre-processed reference data comprises combining the mean    carrier-ambiguity values with interpolated carrier-phase    measurements.-   10. The method of one of 8-9 wherein the acquisition period ranges    from one minute to several minutes and wherein the reference period    ranges from one hour to several hours.-   11. A method of processing data acquired by a GNSS receiver,    comprising:    -   a. Producing a set of pre-processed reference data as in one of        8-10, and    -   b. Preparing single-differenced minimum-error combination        observations between the rover data and the pre-processed        reference data.-   12. The method of 11, further comprising filtering the    single-difference minimum-error combination observations with a    Kalman filter having states for receiver position, receiver clock    and ambiguities.-   13. The method of 11 or 12, wherein the minimum-error combination    observations comprise L-1 phase and an L1-phase plus L1-code    combination-   14. The method of one of 12-13, wherein the Kalman filter has two    rover-clock states and alternates between the two rover-clock states    when processing successive epochs of data.-   15. The method of 14, wherein the Kalman filter has two    rover-position states and alternates between the rover-position    states when processing successive epochs of data.-   16. The method of 14, wherein the Kalman filter has two    rover-position states, further comprising determining from the rover    position states whether the receiver has changed position by more    than a predetermined amount from a prior epoch to a current epoch-   17. The method of 16, further comprising estimating change of rover    position between a prior epoch and a current epoch from phase    measurement differences between the prior epoch and the current    epoch.-   18. The method of 17, wherein estimating change of rover position    comprises applying a decorrelated delta-phase observable.-   19. A method of differential processing of data acquired by a GNSS    receiver with multi-frequency GNSS reference data, comprising    applying a whitening of noise technique to reduce non-Gaussian    positioning errors.-   20. A method of differentially-processing GNSS data in accordance    with one of 1-19, comprising: applying a Kalman filter having    position states for current epoch and for previous epoch and an    ambiguity state for each satellite to prepare a backward estimate of    position and ambiguities for the GNSS receiver, and applying a    Kalman filter having position states for current epoch and for    previous epoch and an ambiguity state for each satellite to prepare    a forward estimate of position and ambiguities for the GNSS    receiver.-   21. The method of 20, wherein the backward estimate and the forward    estimate are prepared with GNSS receiver data collected over a    period substantially longer than a feature-occupation time.-   22. The method of 20 or 21, further comprising combining the    backward estimate and the forward estimate to produced a smoothed    position estimate.-   23. A method of processing data acquired by a GNSS receiver,    comprising:    -   a. Determining a position estimate of the GNSS receiver by        processing the data with reference data from a selected        reference receiver,    -   b. Repeating step a. to obtaining a position estimate of the        GNSS receiver for each of a plurality of reference receivers,    -   c. Averaging the position estimates to produce a final position        estimate of the GNSS receiver.-   24. The method of 23, wherein averaging the position estimates    comprises applying a weighting function proportional to quality of    the position estimates.-   25. The method of 24, wherein quality of the position estimates is    determined from at least one of: number of satellites tracked,    positional dilution of precision, baseline length from the GNSS    receiver to reference receiver, number of epochs of valid data, and    root-mean-square of residual errors.-   26. The method of 24 or 25, wherein quality of the position    estimates is determined from quantities used to estimate covariance    of coordinates of the position estimates.-   27. The method of one of 24-26, wherein the weighting function is in    inverse relation to covariance.-   28. The method of one of 24-27, wherein multiple position estimates    are prepared for each reference station and one of the position    estimates is selected from among the multiple position estimates.-   29. The method of 28, wherein the multiple position estimates    comprise an intermediate L1 float solution and an intermediate L1/L2    float solution.-   30. The method of one of 23-29, wherein the position estimates    comprise intermediate code solutions and wherein averaging comprises    averaging the intermediate code solutions to produce an adjusted    code solution.-   31. The method of one of 23-30, wherein the position estimates    comprise intermediate float solutions and wherein averaging    comprises averaging the intermediate float solutions to produce an    adjusted float solution.-   32. The method of one of 23-29, wherein the position estimates    comprise intermediate code solutions and intermediate float    solutions, wherein averaging comprises averaging the intermediate    code solutions to produce an adjusted code solution and averaging    the intermediate float solutions to produce an adjusted float    solution.-   33. The method of 32, further comprising combining the adjusted    float solution and the adjusted code solution to produce a final    solution.-   34. A method of preparing a post-processing accuracy table,    comprising:    -   a. Acquiring a set of GNSS data,    -   b. Processing a subset of the GNSS data to produce a result for        the subset,    -   c. Repeating step b. for various subsets of the GNSS data to        produce a result for each subset,    -   d. Sorting the results into bins according to common features of        the subset,    -   e. For each bin, determining parameters to characterize the        results of the bin, and storing the parameters and the common        features of the subset in a table.-   35. The method of 34, wherein the parameters characterize results of    the bin as positioning error versus amount of data in the subset.-   36. The method of 34 or 35, wherein the common features comprise one    or more of: number of frequencies of GNSS data, number of epochs of    data, number of satellites tracked, horizontal dilution of    precision, and position dilution of precision.-   37. A method of predicting post-processing accuracy in a GNSS    receiver, comprising:    -   a. Acquiring GNSS data with the receiver,    -   b. Determining at least one feature from the data,    -   c. Using the at least one feature to look up in a        previously-prepared table a set of parameters characterizing a        relationship between positioning error and amount of data        acquired, and    -   d. Determining from the set of parameters at least one of: an        expected post-processed-position error, and an amount of data        needed to achieve a desired post-processed position error.-   38. The method of 37, wherein the parameters characterize results of    the bin as positioning error versus amount of data in the subset.-   39. The method of 37 or 38, wherein the at least one feature    comprises at least one of: number of frequencies of GNSS data,    number of epochs of data, number of satellites tracked, horizontal    dilution of precision, and position dilution of precision.-   40. A stored data structure useful for estimating post-processing    accuracy during collection of GNSS data, comprising: a set of    parameters characterizing a relationship between positioning error    and amount of data acquired for at least one feature of acquired    data.-   41. The data structure of 40, wherein the at least one feature    comprises at least one of: number of frequencies of GNSS data,    number of epochs of data, number of satellites tracked, horizontal    dilution of precision, and position dilution of precision.-   42. A storage medium containing a data structure according to one of    40 or 41.-   43. Apparatus for performing any one or more of the above methods.-   44. Methods and apparatus incorporating a Kalman-filter-based    interpolation scheme that can interpolate reference data to any    desired data rate.-   45. Methods and apparatus apply the L1-only code-plus-carrier    combination for differential positioning.-   46. Super-smoothing of the reference pseudo-range data for    differential positioning (use substantially more reference data than    collected data).-   47. Motion detection using Kalman filter states.-   48. Multi-baseline adjustment (MBA)-   49. Methods and apparatus for predicting post-processing accuracy    while collecting GNSS measurement data, in which parameters of the    collected data are used to look up post-processing accuracy values    in a lookup table prepared from a set of known data.-   50. Combinations of 44-49 to deliver high-precision accuracy with    short convergence time.-   51. Methods and apparatus for predicting post-processing accuracy    while collecting GNSS measurement data, in which parameters of the    collected data are used to look up post-processing accuracy values    in a lookup table prepared from a set of known data.

The invention claimed is:
 1. A method of determining position of a roverfrom observations of GNSS signals, comprising: a. obtaining observationsof GNSS signals at a rover location, b. obtaining observations of theGNSS signals at a plurality of reference stations, each referencestation defining a respective baseline between the rover location and areference station location, c. for each reference station, determining arespective differentially-corrected rover position, wherein at least oneof the differentially-corrected rover positions is based on one of (i) amultiple-frequency ionosphere-free observable combination, (ii) amultiple-frequency code-phase observable combination, and (iii) asingle-frequency carrier-phase and code-plus-carrier-phase combination,and d. preparing a weighted combination of the differentially-correctedrover positions.
 2. The method of claim 1, further comprisingdetermining a quality measure for each of the differentially-correctedrover positions and preparing the weighted combination as a function ofthe quality measures.
 3. The method of claim 2, wherein the qualitymeasure is an inverse covariance.
 4. The method of claim 1, wherein atleast one of the differentially-corrected rover positions is based on asingle-frequency observable combination.
 5. Apparatus for determiningposition of a rover from observations of GNSS signals collected at arover location and at a plurality of reference stations, each referencestation defining a respective baseline between the rover location and areference station location, comprising a. an element to determine, foreach reference station, a respective differentially-corrected roverposition, wherein at least one of the differentially-corrected roverpositions is based on one of (i) a multiple-frequency ionosphere-freeobservable combination, (ii) a multiple-frequency code-phase observablecombination, and (iii) a single-frequency carrier-phase andcode-plus-carrier-phase combination, and b. an element to prepare aweighted combination of the differentially-corrected rover positions. 6.The apparatus of claim 5, further comprising an element to determine aquality measure for each of the differentially-corrected roverpositions, wherein the weighted combination is prepared as a function ofthe quality measures.
 7. The apparatus of claim 5, wherein the qualitymeasure is an inverse covariance.
 8. The apparatus of claim 5, whereinat least one of the differentially-corrected rover positions is based ona single-frequency observable combination.
 9. An article of manufacturecomprising a computer-readable medium having instructions stored thereinfor causing a machine to perform the method of claim 1.